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^3 ! Abstract. The thermal conductivity of doped graphene flake of finite size is 

I ■ investigated with emphasis on the influence of mass of substituting atoms on this 

property. It is shown that the graphene doping by small concentrations of relatively 
heavy atoms results in a disproportionately impressive drop of lattice thermal 
conductivity. 
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1. Introduction 

This paper was motivated by recent investigations p!JI21|3j|IJ[5j|6j[71[8j[9], * n wn i° n a 
real possibility of the experimental measurement of the thermal conductivity of graphene 
was established and its extremely large value was confirmed. 

The proposed and performed experiments were based there on the effect of the 
temperature shift of the Raman G peak position due to the anharmonicity of graphene 
lattice pm E]. Since that shift may be measured with sufficient accuracy, the thermal 
conductivity of laser heated graphene may be restored using the heat equation with 
known sources and given geometrical parameters of the graphene flake. Bearing in mind 
the efficacy of existing experimental methods we consider in this paper the influence of 
moderate concentrations of different point defects (carbon isotopes, substituting atoms 
and chemically adsorbed molecules) on the thermal conductivity of graphene. 

Note that real and numerical experiments showed that graphene can be doped if 
required with boron, nitrogen (T2J [13], aluminium [TJ] etc with concentrations reaching 
as much as several per cent without breaking the graphene planar structure [151 US] • 

The main task of this work is the elucidation of how the thermal conductivity 
of doped graphene depends on the temperature and mass of defects using as a tool 
classical methods of the quantum theory of solids. In doing so we restrict ourselves by 
modeling all the substituting atoms as isotopic defects, assuming that the account of 
bond weakening may only give an additional reduction of lattice thermal conductivity 

This paper is organized as follows. In Sections 2 we discuss the main mechanisms 
of the thermal resistance in graphene. In Sections 3 and 4 we describe briefly the effect 
of phonon-phonon scattering and edge scattering as the dominant thermal resistance 
mechanisms in the case of ideal graphene. In Section 5 the explicit expression for the 
phonon relaxation time determined by phonon-defect scattering is deduced as a function 
of the quasi-momentum for each phonon branch. In the last section we summarize 
the obtained results in the form of a fitting formula for the thermal conductivity of a 
graphene flake as function of temperature, flake size, mass and concentration of dopant. 

2. Heat transport in graphene and mechanisms of thermal resistance 

The heat carriers in graphene are the phonons and conductive electrons. The ideal 
graphene is a semi-metal with massless charge carriers. As follows, electron and phonon 
contributions to the heat transport should be in the ratio ~ c p h./c e i where c p h. and 
c e i are characteristic velocities of acoustic phonons and electrons (holes) in graphene. 
Due to the near 50-fold difference in velocities of electrons and phonons, the electron 
contribution to the ideal graphene total thermal conductivity appears to be less than 
1% |2j. Besides in semi-metals the number of carriers at low and room temperatures is 
small compared to that of lattice atoms. Since the phonon scattering on free carriers 
is proportional to their concentration, then the effect of phonon-electron scattering on 
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the lattice heat transport is small comparing to that of phonon-phonon and phonon- 
defect scattering |T7] . Therefore further in this paper we will assume that the electronic 
subsystem counts very little in the heat transfer as a whole and can be ignored, the more 
so that the error of present measurements of graphene heat conductivity is > 10%. With 
this assumption the finite value of the thermal conductivity of real graphene results from 
phonon scattering on other phonons, flake edges and point defects. 

Phonon-phonon scattering is conditioned only by the anharmonic contributions 
to the potential energy of inter-atomic interaction. The quasi-momentum and energy 
conservation laws only allow phonon scattering processes in which at least three phonons 
are involved. As the probability of collisions of four or more phonons is rather small, 
especially at low temperatures where the equilibrium number of phonons is not great, 
their contribution to the thermal resistance may be ignored. Remind that multi-phonon 
scattering processes are divided into normal processes (N-processes) with the strict 
conservation of momentum and so-called umklapp processes (U- processes), for which 
the conservation of total quasi-momentum holds up to a non-zero reciprocal lattice 
vector. The normal processes do not change the total heat flux, thus they drop out 
from possible mechanisms of thermal resistance, while the U-processes contribute, in an 
important manner, to the thermal resistance of semiconductors and isolators at high 
temperatures [T7j. On the other hand, with decreasing temperature the contribution 
of U-processes to the thermal resistance of graphene is suppressed by effects of the 
defect-induced and edge phonon scattering. 

Considering phonon-defect scattering, we restrict ourselves to study only the impact 
of isotopic defects on the heat transfer in graphene. Note that along with 13 C isotopes 
(present in natural carbon at a concentration of ~ 1%) the chemically adsorbed atoms 
and substitutional atoms (like boron, nitrogen, aluminium etc) can be to some extent 
treated as isotopic defects. In reality, concentrations of all mentioned defects may range 
from zero up to several per cent. However, as we show below, for the low atomic mass of 
carbon even a concentration of heavy "isotopic" defects < 1% may double the thermal 
resistance. 

At low temperatures the intensities of all mentioned processes rapidly decrease and 
the edge scattering is thus rendered the main mechanism limiting the heat transport. 
Taking account of the very high Debye temperature of graphene and exponential drop 
in the probability of U-processes there is a temperature range within which only edge 
and defect scattering should be considered. 

Taking into account all the above mentioned contributions, the total relaxation 
time of a certain phonon state is given by 
1111 
T Ks r Ks r k J r k J 

where k is the phonon wave vector and s is the oscillation branch; r^ s , t*/ and r^ s 9e 
are the relaxation times determined by the phonon-phonon, phonon-defect and edge 
scattering, respectively. 
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If the temperature of graphen is well below its Debye temperature (the latter is 
near 2000K), then the phonon-phonon interaction can be considered as a rather small 
perturbation and thus different phonon states can be considered as independent. Thus 
the total lattice thermal conductivity k can be represented as a sum over all possible 
phonon modes [17J 



K = g C k,^k,sAk,s = g C M%,s r k, s , (2) 



.2 

k,s k,s 

where k runs the first Brillouin zone, Ck, s is the specific heat per normal mode, 

' hu s (k)\ 2 1 
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where s numbers the phonon branches, u^ s is the frequency of the |k, s)-phonon, 
Vk, s = Vk^k.s is the phonon group velocity and Tk )S is defined by (JT]). 

Since our main task is to study the impact of defects on the thermal conductivity of 
ideal graphene, we will not scrutinize the contribution of the phonon-phonon U-processes 
using the same general expression as in (2j HJ [TBI US]; however, with w a (k) defined as 
described below <^}: 

Tk - (T) = ^T^T^k)' S = 1 ' 2 ' 3 ' ( 3 ) 

where m is the atomic mass of 12 C isotope and ujd jS are conventional Debye frequencies 
for three acoustic phonon modes (wd,i — 2.66 x 10 14 , ujb,2 = 2.38 x 10 14 and 
w Di 3 = 1.32 x 10 14 rad/s). 

In order to find from ([2]) certain values of k we use the explicit expressions for 
Uk, s in graphene obtained in (SS] on the basis of the following computationally simple 
harmonic nearest neighbor model. 

Considering the total potential energy W of ideal graphene as a sum of three 
components depending on different types of interaction between the nearest neighbors: 

- in-plane central forces with a force constant J\\ 

- in-plane non-central forces with a force constant J2; 

- out-of-plane non-central forces conditioned by 7r-electrons with a force constant 

Js, 

and applying Bloch's theorem one can reduce within the framework of this model the 
problem of finding the sought frequencies u s (k) to evaluation of eigenvalues of the ideal 
graphene dynamical matrix 



where 




D<k) = ^{WD*) < 4) 
/ Ji + 2J 2 \ 



Ji + 2J 2 

\ 2j 3 / 
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In this way we find three acoustic and three optical branches 0J s (k) of phonon spectra 
(figured]). With notation 

Xq = J1J2 + Xi = Jf + I6X0, X2 = Jf — 8X0, 

F (k) = 2 [cos {kia — k,2a) + cos k\a + cos/^a] , 

F x (k) = 12(J 2 + 2X ) + ^(k)^ 2 + 8X ), (5) 
F 2 (k) = {9 J 2 + Xi + 2X 2 cos (k 1 a — k 2 a) (cos k\a + cos k 2 a) + 

+2 cos kia cos k 2 a [X 2 +Xi cos {kia — k 2 a)] — 3J 2 -F (k)} 5 , 

the expressions for the eigenfrequencies w s (k) take form: 

1/2 

(6) 



uzA,zo(k) 



3± J3 + F (k) 



VLA,TA,LO,ToM = ± il W ± 2 (k) 



1/2 



(7) 



The values of the force constants 
Ji = 135 J in" 2 , 
J 2 = 245 Jm" 2 , 
J 3 = 83.9 Jm" 2 . 

were chosen in [20] to fit in the experimental data on inelastic X-ray scattering in 
graphene [21]. It is noteworthy that the model with these Ji, J 2 , J 3 reproduces Raman 
experiment results and also gives appropriate values for sound velocities 

vla = 18.4 km s -1 , vta = 16.5 km s _1 , vza= 9.2 km s -1 (8) 

in graphene (see [2T] 122] 123] ) . 

We remark that the given expression (]6]) for the out-of-plane (bending) mode uza(X) 
is inconsistent in the limit |k| — > with known dependence uzaO*-) ~ k 2 . However the 
quasi-linear dependence of uzaO<.) — vza |k| which follows from ([6]) occurs in graphene 
w^ y 4(k)-dispersion curves obtained in different ways already for |k|>^r-^--^Tor for 
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^za ^2-^4 meV [U [21] |2U [25]. This gives grounds for using the formula ([6]) for 
T > 20 — 30 K. It is worth recalling that the out-of-plane phonon mode add little to the 
thermal transport in pure graphene because of its low group velocity and large Gruneisen 
parameter [3] H] [9] [26]. By our calculations, the contribution of this mode to the total 
graphene thermal conductivity at room temperature does not exceed ~ 15 — 20%. 

In all our calculations we used the exact expressions (j6]) without any further 
simplification. At low temperatures (< 100K) this does not improve visibly the Debye 
approximation. However, the discrepancy between results of the Debye approximation 
and those obtained using ([6]) are quite visible even at room temperature, especially for 
small graphene flakes (figure [2]) . 
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3. Effect of the edges 

In the low temperature limit where the U-scattering processes are not active and defect 
scattering is relatively weak the phonon mean free path A grows rapidly with T~ l . On 
the other hand A cannot exceed the lattice linear size. Thus the thermal conductivity 
in this region is limited only by the edge scattering and should depend linearly on the 
scaling factor for similar graphene flakes. 

For a bounded convex graphene flake with diameter much greater than the lattice 
constant, the distribution of eigenfrequencies is just the same as that in the case of 
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Figure 2. Thermal conductivities of ideal graphene flakes of different sizes 
calculated using the exact expressions (j6|) (solid lines) and obtained within the Debye 
approximation (dashed lines). 
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infinite ideal graphene. The phonon eigenfrequencies of the flake numbered in increasing 
order are well within the dispersion curves of ideal graphene. Though the displacement 
eigenvector |k, s) for the infinite ideal graphene lattice is not an eigenvector for a 
bounded graphene flake, it is, however, the superposition of the flake's displacement 
eigenvector corresponding to frequencies close to Uk yS . Therefore the uncertainty of 
the exact eigenfrequencies can be reduced to the uncertainty of quasi-momentum k: 
Alu = ffAA;, where AA; ~ = S is the flake's area and L = \/S is its characteristic 
size. 

For the phonon wave associated with |k, s) the energy-time uncertainty principle 

r £ £AE > f g ives 



Setting 



< AE = h—Ak ~ h— — = 2ixh^. 



2r^ e ~ dk dk L 



2 IT 



A ed9e = ^sj dxd yj v&y^W = ( 9 ) 

S 

where D(x,y,4>) is the length of the segment inside the flake of ray coming from the 
point (x,y) and ^-inclined toward the x-axis, and / is the dimensionless flake form-factor 
(/ = 5^§72~ 0.479 for the disc and reduces with increasing of the flake diameter) we use 
further for the r^ e an approximate formula 

/ edge\-l _ v k,s _ v k,s / ln x 
> - A edge ~ ^1 

in which the ballistic mean free path A edge absorbs all the peculiarities of the flake 
geometry. 

In the case of graphene with abundant linear structural defects the boundaries of 
emerging domains became the main source of the edge scattering. If so, to obtain L one 
may calculate it as above for a certain domain with some shape and area Sd and then 
average the obtained value over the distribution of the domain shapes and areas and 
also incorporating in this case the specularity parameters p at the domain boundaries 
as in [3]. 

Note that at low temperatures the finite size of a graphene flake is a dominant 
factor dictating the finiteness of the phonon free path. Since all the active phonons are 
acoustic phonons near the T-point where u s ~ v s k, then the thermal conductivity takes 
the form k = \k ed9e v s C s ~ T 2 , where s numbers acoustic phonon branches and C s 

s 

are corresponding partial lattice specific heats per unit volume. 
4. Defect scattering 

The mean free path of the |k, s)-phonon caused only by its scattering on point defects 
is defined by the evident expression 

\def = 1 

k ' s n<r 8 (k)' 
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where n is the surface density of point defects and cr s (k) is the corresponding phonon 
total cross-section. Hence 
1 



def 



na s (k) |v Sjk |. (li; 



Using the optical theorem and Lippmann-Schwinger formula for calculation of <r s (k) 
in dTTJ) yields [27] (see Chapter 11) 

1 il Im(k,s|T+GT|k,s) = — Im(k,s|T|k,s). (12) 



where Q is the area of a two-dimensional graphene lattice, 

G = G(k, u 2 + i 0) = (D(k) - I(u 2 + iO)y X 

is the Green's function of the infinite ideal graphene lattice with the dynamical matrix 
D(k) defined above (jl"]) , and T = T(w£ s -MO) is the T-matrix, which connects the 
Green's function G of a defect lattice with that of the ideal lattice: G = G — GTG, or 
T = D(l + GD) 1 . For an isotopic defect located at the lattice site (0,0) the element 
of T-matrix with indices 7, j, I"; 7', j', I'" (where 7 numbers graphene sublattices, j is 
associated with spatial components of the atom displacement and I is the lattice vectors) 
has the form 

TifUl", V") = -w 2 ^ 7 ,o *i« (Q- 1 )"" S Yt0 8 v „ i0 , 



and 



lm(u,k\T(ul + iO)\u, k) 



(13) 



Im (col + * 0) E «„(*:))' ([Q. 2+i0 ] _1 ) M e^(k). 



1 Am 

N m 

hJ 



(14) 



The expressions (fT2|) . (jTBl . (|T4|) were immediately used for the calculation of 



5. Discussion 

We used directly the relations for the total lattice thermal conductivity (II]), ((2]) together 
with the explicit expressions discussed above ([3]), ({TO]) . ( TT2]) and (J6]) for r^ s , t^^ 6 , r^ 6 / 
and tUk )S , to calculate the thermal conductivity of graphene flakes in a wide range of 
linear sizes, temperatures and concentrations of isotopic defects with n = — varying 
from to 2. 

The obtained data on graphene thermal conductivity in the range T e (20, 500)K, 
L G (1, 100) /im, < 2 and n < 0.05 are reproduced with an accuracy to < 1% by 
the following empirical formula which permits us to obtain the value of the thermal 
conductivity without re-calculating explicit expressions ([3]), ffTO]) . f|T2]) and (any 
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further improvement of this formula has no sense in view of the approximations laid 
down into the above theoretical constructions) 

-l 



k(Q, C, /i, n) 



1 



fi 2 n 



(15) 



kj(&,C) K D (Q,C,n,n 

where Ki(Q,C) is the thermal conductivity of an ideal graphene flake, G = T/6 
6 = hcuD^kg 1 and C = is the characteristic length of the flake in micrometers, 



and 



A 1 (C) = 



[6 + A 2 (C)f + A&C)' 
1310 . 3 



0.573 + £-o-45 



26 C 



— - 0.0594; A 3 (C) 



121 £ 



0.35 



- 0.005; 



k d (Q, C, fi, n) = Bx{C) V6 + B 2 (C)fM n°- 5d + B 3 (C) 
Ih{C) = -14.9 £ 7 / 13 + - - - - 107, B 2 (C) = 46.5^^ + 913, 



11 - log £ 
B 3 {£) = 13.6 C 0A2 - 3.7 log C - 16. 
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Figure 3. Thermal conductivity for square graphene flakes of different characteristic 
sizes with aluminium-27 and nitrogen-14 (or carbon-14) isotopic defects. 



As expected, the results for ideal graphene are close to those of Balandin et al. [1] . 
Notice that the qualitative trend of the thermal conductivity of doped graphene as a 
function of temperature, type and concentration of defects emerging from (I15j) agrees 
well with arguments and qualitative conclusions for graphene of the comprehensive 
review paper [28]. At low temperatures the dominant factor of phonon scattering is 
the finiteness of a graphene flake that gives k ~ T 2 . For relatively high temperatures 
the thermal conductivity k varies in inverse proportion to T due to the phonon-phonon 
scattering. The phonon-phonon scattering is calculated here by summation over all the 
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Brillouin zone in contrast to [31 S] , where the Debye approximation was used. The Debye 
approximation tends to increase the value of thermal conductivity since it gives non-zero 
sound velocities at the Brillouin zone edges. The difference between these approaches is 
noticeable for T > 100 — 200K being as much as 20% for T = 500K, but so far remaining 
within the experimental error. 

In order to avoid any additional geometrical and non-geometric parameters, the 
phonon mean free path associated with the edge scattering is taken simply as fl9]). 

The established decrease of graphene thermal conductivity caused by isotopic 
defects appears to be approximately linear in concentration n (for small concentrations) 
and strongly dependent on fi and temperature. Note that though the cross-section 
of phonon scattering on defects is by itself temperature independent, it depends 
substantially on the phonon quasi-momentum. Therefore the relative contributions of 
phonons from different parts of the first Brillouin zone to the heat transport changes 
with temperature, resulting in the complex temperature dependence of the contribution 
of defects to the total heat resistance. 

For substitution defects with \fi\ < 1/6, such as boron, nitrogen and carbon isotopes 
the relative reduction in thermal conductivity is just proportional to n/x 2 and is well 
within experimental error. For example, for thermal conductivity of graphene with 
concentration n ~ 0.01 of 13 C according to the above results at room temperature the 
derivative |^ ~— 1 (WmK -1 ) and varies from -1.2 to -0.8 as L increases from 1 to 100 
jum. For the isotopically modified 13 C-graphene with concentration n ~ 0.01 of 12 C this 
derivative is just the same while the difference in heat conductivities of pure 12 C- and 
13 C-graphene appears to be lesser than that between pure 12 C- and pristine graphene. 
This agrees with experimental results on the thermal conductivity of graphene with 
concentrations of 13 C isotope of 0.011 and 0.992 [TJ. Furthermore, despite the fact 
that the model developed above is not valid without further improvement for very high 
defect concentrations, the ratio of the calculated value of thermal conductivity of pristine 
graphene to that of 50-50 isotope mixture is the same as in [7J. In the meantime we 
cannot give so far the satisfactory explanation of the drastic drop in the measured value 
of graphene thermal conductivity for the ~1% isotope concentrations observed in [7J. 

For heavier defects (aluminium atoms and different chemically adsorbed molecules) 
with // > 1 the thermal conductivity drops sufficiently even for n ~ 1%. As an example, 
1% of "isotopic" 27 A1 (/i = 1.25) reduces the thermal conductivity by more than half 
(figure [3]). Strictly speaking, the replacement atoms and adsorbed molecules are not 
isotopic defects. However, the entailed bonds weakening may only enhance the effect 
of low frequency phonon scattering (see [17J, Chapter VI) and thus lead to a further 
reduction in the thermal conductivity of graphene. 

It should be noted that chemisorbed molecules may desorb in the temperature 
range 300-600 K decreasing the concentration of defects. So the graphene thermal 
conductivity under gas-surface equilibrium conditions in corresponding temperature 
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regions may show a sufficient rise with T instead of the "normal" behavior k ~ T _1 0. 

Thus the insertion of rather low concentrations of rather heavy dopants has a 
disproportionately great impact on the thermal resistance of graphene. This conclusion 
is completely in line with results of rather recent molecular dynamics simulations of 
phonon transport in doped graphene [281 1291 30]. 
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